Modal Analysis of a AA class Violoncello spruce top plate

This live script presents the complete workflow followed by the authors in the process of determining the dynamic properties of a Violoncello top plate.
After making use of Matlab's own example on the modal analysis of a flexible wing aircraft for the development of this code [1], the authors felt it necessary to share this example in the hopes that it would help others seeking to perform experimental modal analysis through Matlab. Said example was born as a result of one of the authors master's thesis which focused on the manufacture of classical instruments using advanced materials. And was made possible thanks to the help of spanish luthier Carlos Moreno (www.luthieropera.com) who lended the detached top plate used during tests.
As a consequence, the dynamic properties of the currently accepted as a "good" instrument were measured and set as the goal to reach using advanced design and manufacturing techniques. The modal assurance criterion is seen as a way of comparison between not only the different test results but also with test from a simulation.
Another extra feature is the capability of this code of plotting the mode shape using the entire set of 3D mesh points, although only 16 of them had been measured. This is achieved through polynomial fitting on the mode shape components and is application to the rest of the plate's geometry. The authors wish to note that the method used here is only valid for 2D geometries and think that a nice future work would be to extend this to any geometry.

The Violoncello

Parts

Although the Violoncello is comprised of a series of components, all of them involved in the generation of sound (See "the physics of the violin" by Lothar Cremer for further reference), the majority of the sound pressure waves emanate from the top plate-f holes. Therefore, in order to recreate the instrument's characteristic timbre using non-conventional materials dynamic properties of the top plate should be replicated as a minimum.

Obtaining 3D geometry

A 3D scan of both sides of the top plate was performed using an affordable secondhand scanning device. After some treatment of both surfaces, a merger between them was possible creating a 3D geometry.
Creating a Mesh
A mesh was then associated to the 3D geometry using freeware and imposing symmetry between the left and right halves of the geometry. The mesh associated with the longitudinal reinforcement was created separately and restrictions on displacements between them both were imposed.
Test DOFs location
Using other parametric CAD software an array of 16 test points was selected ensuring their symmetric placement.
A template for DOF location on the top plate was made by printing a 2d projection of said mesh.

Test setup

Accelerometers placing

Using the wax provided by the accelerometer manufacturer, accelerometers were attached to the top plate and their DOF numbers were marked on the cables.

Hardware / software

A NI 9185 chassis was used together with NI 9234 data acquisition cards for both the accelerometers and the impact hammer signals. The selected accelerometers had a nominal sensitivity of 100mV/g (although actual calibrated sensitivities were input into the NImax task). and the impact hammer had the lowest possible mass (250gr) and a rubber tip in order to minimize damage to the plate's surface. All test files were stored as ascii formatted data with headers and were converted into .mat format for the sake of simplicity of this live script.
The instrumented plate was hung from a structure using very flexible rubber bands in order to ensure close to free-free boundary conditions during the tests. And the posterior parts of the accelerometers were marked with tape in order to impact on the positive sense of each accelerometer by hitting the plate from its rear side.

Test sequence

The test sequence is summarized in the following table.

Matlab code

Importing test results into Matlab

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Violoncello Experimental Modal Analysis sequence %%%
%%% Authors: Álvaro Nieto | Paula González %%%
%%% Date: January 10th 2020 %%%
%%% Converted to live script on April 9th 2020 %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clearvars
%% Loading Data
load Modal_tapa_chelo_abeto_HAMMERS_Y_SENSORS.mat
% Common parameters setting
fs = 10240; % data sampling frequency (Hz)
Ts = 1/fs; % sample time
%
%Creating list of hammer signals
hammer_list=whos('hammer_*');
sensors_list=whos('sensors_*');

Setting accelerometer positions and orientation and loading 3D model geometry for improved mode shape visualization

%% Accelerometer positions & orientations (reference coordinates frame)
AccelPos = [...
0.000 0.745 0; %DOF 1
-0.103 0.618 0; %DOF 2
0.103 0.618 0; %DOF 3
0.000 0.557 0; %DOF 4
0.000 0.454 0; %DOF 5
0.050 0.320 0; %DOF 6
0.000 0.270 0; %DOF 7
-0.108 0.252 0; %DOF 8
0.108 0.252 0; %DOF 9
0.000 0.212 0; %DOF 10
-0.137 0.168 0; %DOF 11
-0.087 0.168 0; %DOF 12
0.087 0.168 0; %DOF 13
0.137 0.168 0; %DOF 14
-0.092 0.088 0; %DOF 15
0.092 0.088 0]; %DOF 16
AccelDir = [...
0 0 1; %DOF 1
0 0 1; %DOF 2
0 0 1; %DOF 3
0 0 1; %DOF 4
0 0 1; %DOF 5
0 0 1; %DOF 6
0 0 1; %DOF 7
0 0 1; %DOF 8
0 0 1; %DOF 9
0 0 1; %DOF 10
0 0 1; %DOF 11
0 0 1; %DOF 12
0 0 1; %DOF 13
0 0 1; %DOF 14
0 0 1; %DOF 15
0 0 1]; %DOF 16
% Loading 3D points for improved animation
load('Violoncello_Top_plate_XYZ.mat');
%

User selected test visualization, signals vs timestamp

Note both the saturation of the accelerometer signals at 50g during those tests where the excitation force was excessive and the partial detachment of some accelerometers after impact.
Testnumber=11;
post=0.3;%time width display in seconds
switch Testnumber
case num2cell(1:4);DOF=6;case num2cell(5:8);DOF=5;case num2cell(9:12);DOF=14;case num2cell(13:16);DOF=3;
end
Hammer_EMA_test_visualization(Testnumber,eval(sensors_list(Testnumber).name),'g',eval(hammer_list(Testnumber).name),'N',DOF,Ts,post)
hold
Current plot held

Creating iddata objects

The dynamic properties of the plate are modeled by the frequency response functions of the measured degrees of freedom. Said FRFs are identified later on using Matlab's etfe function, which makes use of the test data packed into iddata objects. The following code merges results from impacts on the same DOF into single iddata objects.
%% Creating individual test iddata
%DOF #6
iddata_dof6=merge(iddata(sensors_01,hammer_01,Ts),iddata(sensors_03,hammer_03,Ts)); %Only valid "clean" data included
%iddata_dof6=merge(iddata(sensors_01,hammer_01,Ts),iddata(sensors_02,hammer_02,Ts),iddata(sensors_03,hammer_03,Ts),iddata(sensors_04,hammer_04,Ts));
%DOF #5
iddata_dof5=merge(iddata(sensors_05,hammer_05,Ts),iddata(sensors_06,hammer_06,Ts),iddata(sensors_07,hammer_07,Ts),iddata(sensors_08,hammer_08,Ts));
%DOF #14
%iddata_dof14=merge(iddata(sensors_09,hammer_09,Ts),iddata(sensors_10,hammer_10,Ts),iddata(sensors_11,hammer_11,Ts),ddata(sensors_12,hammer_12,Ts));
iddata_dof14=merge(iddata(sensors_10,hammer_10,Ts),iddata(sensors_11,hammer_11,Ts),iddata(sensors_12,hammer_12,Ts));
%DOF #3
iddata_dof3=merge(iddata(sensors_13,hammer_13,Ts),iddata(sensors_14,hammer_14,Ts),iddata(sensors_15,hammer_15,Ts),iddata(sensors_16,hammer_16,Ts));
%

FRF estimation

A total of 8192 data points is chosen by default as it is a power of two and Matlab will make use of the FFT algorithms. Furthermore, at a data acquisition rate of 10240 Hz, this number of data points represent close to 0.8 seconds of time data, more than enough to capture the entire free decay vibration of all accelerometers (see the previous plots).
%% FRF estimation using etfe Nfpoints Hanning window data points
Nfpoints=8192; % data points for FRF estimation
%Impacts on DOF#6
FRF_DOF6=etfe(iddata_dof6,[],Nfpoints);
FRF_DOF6.OutputName='Acc';%Acc(1),Acc(2),...,Acc(16)
FRF_DOF6.InputName='DOF 6';
FRF_DOF6.InputUnit='N';
FRF_DOF6.OutputUnit=['g';'g';'g';'g';'g';'g';'g';'g';'g';'g';'g';'g';'g';'g';'g';'g'];
%Impacts on DOF#5
FRF_DOF5=etfe(iddata_dof5,[],Nfpoints);
FRF_DOF5.OutputName='Acc';%Acc(1),Acc(2),...,Acc(16)
FRF_DOF5.InputName='DOF 5';
FRF_DOF5.InputUnit='N';
FRF_DOF5.OutputUnit=['g';'g';'g';'g';'g';'g';'g';'g';'g';'g';'g';'g';'g';'g';'g';'g'];
%Impacts on DOF#14
FRF_DOF14=etfe(iddata_dof14,[],Nfpoints);
FRF_DOF14.OutputName='Acc';%Acc(1),Acc(2),...,Acc(16)
FRF_DOF14.InputName='DOF 14';
FRF_DOF14.InputUnit='N';
FRF_DOF14.OutputUnit=['g';'g';'g';'g';'g';'g';'g';'g';'g';'g';'g';'g';'g';'g';'g';'g'];
%Impacts on DOF#3
FRF_DOF3=etfe(iddata_dof3,[],Nfpoints);
FRF_DOF3.OutputName='Acc';%Acc(1),Acc(2),...,Acc(16)
FRF_DOF3.InputName='DOF 3';
FRF_DOF3.InputUnit='N';
FRF_DOF3.OutputUnit=['g';'g';'g';'g';'g';'g';'g';'g';'g';'g';'g';'g';'g';'g';'g';'g'];
%

Estimated FRF plots

In order to get an idea of the number of modes present in the top plate, a Bode plot of some of the DOFs is useful in order to assess both the number of modes present in the system and any areas of the FRFs that might denote poor estimation.
DOF=3; %Plotting the FRF estimated by hitting in one of the tested DOF
Output_DOFs=[1 3 6]; % Select a series of FRF_ij to be plotted, if left empty all available FRFs will be plotted
freq_range=[0 5000]; % Frequency range in Hz for Bode plots of the FRFs, if left empty the entire range will be plotted
plot_FRF_Bode(eval(strcat('FRF_DOF',num2str(DOF))),Output_DOFs,freq_range)

Checking test data quality (coherence function)

After a brief view of the previous plots, one can easily identify at least 20 resonance candidates. However, since the impacts were made using a rubber tip and the accelerometers were fixed using wax, one should be suspicious of the high frequency contents in the identified FRFs. One good way of assessing the test quality is the coherence function [2], which indicates linearity between the frequency components of the input signal and the output signal.
%%
clf
Testnumber=7;
Output_DOF=5;
freq_range=[0 5000];
switch Testnumber
case num2cell(1:4);DOF=6;case num2cell(5:8);DOF=5;case num2cell(9:12);DOF=14;case num2cell(13:16);DOF=3;
end
[frf,f,coh] = modalfrf(eval(hammer_list(Testnumber).name),eval(sensors_list(Testnumber).name),fs,Nfpoints);
fig=figure('Name',strcat('FRF quality Test#',num2str(Testnumber)));
set(fig, 'Position', [0 0 1000 500])
yyaxis left
plot(f,10*log10(abs(frf(:,Output_DOF))))
if isempty(freq_range)
else
xlim(freq_range)
end
grid on
grid minor
yyaxis right
plot(f,coh(:,1))
if isempty(freq_range)
else
xlim(freq_range)
end
title(strcat('FRF (dB) & coherence test#',num2str(Testnumber),' Output DOF #',num2str(Output_DOF),' excitation on DOF#',num2str(DOF)))
The plot shows good coherence (close to 1) between the output DOF acceleration and the input hammer force for FRF components below 3-3.5kHz (depending on the plotted DOF). Although better results could be achieved using adhesives or screws to fix the accelerometers and metal tips on the impact hammers, these options would almost surely damage the plates surface or even crack the plate itself.
Further more, the pure tonal range of the instrument goes from a C2 (65.41 Hz) to A5 (880 Hz) with the highest harmonic achievable by a skilled musician being A7 (3520 Hz). Hence it would be relatively safe to say that the obtained data can be used for representation of the majority of the frequency content responsible of the pure tone response of the plate. Please note that when one note is played, the resonance within the instrument causes a series of harmonic components that couple with the played note amplifying the sound pressure. The amount and amplitude of the harmonics is what musicians usually refer to as the timbre or tone color of the instrument
Those harmonics can have frequencies that go well beyond the played frequency, so in order to recreate the sound emitted by the plate with maximum fidelity its entire FRF frequency range should be replicated. Obviously, this has some practical implications that make the task at hand easier said than done.

Dynamic system estimation, using n4sid or ssest

This section deals with the dynamic system identification from the test data located in the iddata objects. The default state-space model fits the FRFs quite well in the 20–1000 Hz region, hence including the playable C2-A5 range previously explained. We used a 50th-order model which means that a maximum of 25 oscillatory modes can be identified with the fitted system.
The users are encouraged to extend/reduce the frequency range and the state space order in order to get a better understanding of the calculation times required and the quality of the achieved results.
Furthermore, an iterative algorithm can be used making use of the dropdown menu, the authors might like to point out that although this method greatly improves results with a large number of iterations, it comes with the cost of long computation times.
DOF=5; %Estimation on the iddata object corresponding to the selected DOF hits tests
SS_Order=50; %State space order, note that number of system modes will be a maximum of SS_Order/2
iterative_search='yes' ; %non iterative linear solver if 'no', Levenberg-Marquardt search method if 'yes'
Max_Iterations=3; % Maximum number of iterations performed by the iterative system estimation fit
freq_range=[20 1000]; %Frequency range of the FRF functions to be fitted
f = eval(strcat('FRF_DOF',num2str(DOF),'.Frequency/2/pi')); % extract frequency vector in Hz (G stores frequency in rad/s)
Gs = fselect(eval(strcat('FRF_DOF',num2str(DOF))),freq_range(1)*(2*pi),freq_range(2)*(2*pi)); % "fselect" selects the FRF in the requested range
FRF = squeeze(Gs.ResponseData); % Takes FRF array from the selected range
Weighting = mean(1./sqrt(abs(FRF))).'; % Weighting applied as inverse of the squared FRF value, so that
switch iterative_search case 'no';iterative=0;case 'yes';iterative=1;end
% low value FRF points such as antiresonances are modeled correctly
if iterative==0
n4Opt = n4sidOptions('EstimateCovariance',false,...
'WeightingFilter',Weighting,...
'OutputWeight', eye(size(eval(strcat('FRF_DOF',num2str(DOF))),1)));
sys_DOF_Hits = n4sid(Gs,SS_Order,'Ts',0,'Feedthrough',true,n4Opt);
elseif iterative==1
n4Opt = n4sidOptions('EstimateCovariance',false,...
'WeightingFilter',Weighting,...
'OutputWeight', eye(size(eval(strcat('FRF_DOF',num2str(DOF))),1)));
ssOpt = ssestOptions('EstimateCovariance',true,...
'WeightingFilter',Weighting,...
'SearchMethod','lm',... % use Levenberg-Marquardt search method
'Display','on',...
'OutputWeight',eye(size(eval(strcat('FRF_DOF',num2str(DOF))),1)));
ssOpt.SearchOptions.MaxIterations = Max_Iterations;
sys_DOF_Hits = n4sid(Gs,SS_Order,'Ts',0,'Feedthrough',true,n4Opt);
sys_DOF_Hits = ssest(Gs, sys_DOF_Hits, ssOpt); % estimate state-space model (this takes several minutes)
end
Warning: Information content in input data does not support estimation of initial states. The "X0" property will be set to zero.
Fit = sys_DOF_Hits.Report.Fit.FitPercent';
[~,Imin] = min(Fit); %The worst fit
[~,Imax] = max(Fit); %the best fit
%Plotting "best" and "worst" results
opt = bodeoptions; % plot options
opt.FreqUnits = 'Hz'; % show frequency in Hz
opt.PhaseVisible = 'off'; % do not show phase
fig=figure
fig =
Figure (LiveEditorRunTimeFigure) with properties: Number: 4 Name: '' Color: [0.9400 0.9400 0.9400] Position: [488 342 560 420] Units: 'pixels' Show all properties
set(fig, 'Position', [0 0 1000 500])
bodeplot(Gs([Imin, Imax],:), sys_DOF_Hits([Imin, Imax],:), opt);
xlim(freq_range)
if iterative==0
title(strcat('Worst (top) and best (bottom) fits between data and model, hits on DOF ',num2str(DOF)))
elseif iterative==1
title(strcat('Worst (top) and best (bottom) fits between data and refined model, hits on DOF ',num2str(DOF)))
end
grid on
grid minor
legend('Data', 'Model')

Modal parameters from fitted model

Once a dynamic system is identified, the following parameters are extracted using matlab's modalfit function:
%% Modal parameter computation, using fitted dynamic system from hammer hits
%Use the modalfit command to fetch the natural frequencies in Hz for these modes.
[fn_sys_DOF_Hits,dr_sys_DOF_Hits,ms_sys_DOF_Hits] = modalfit(sys_DOF_Hits, f, SS_Order/2); % natural frequencies (Hz)
Warning: 'MaxModes' value exceeds the maximum modal order. Defaulting to 24 modes.

Modal animation

This section plots the selected mode shape along with its natural frequency and damping ratio. Select a mode number and press CTRL+ENTER keys to animate it.
mode=1; %Animate by selecting the desired mode and pressing CTRL+ENTER keys
AnimateOneMode(mode,fn_sys_DOF_Hits,dr_sys_DOF_Hits,ms_sys_DOF_Hits,AccelPos,AccelDir,VioloncellotopplateCADcoordinates)
%
Modal Assurance Criterion
This section computes and plots the modal assurance criterion for the identified dynamic system modes (auto-MAC). Ideally non coupled modes properly identified should produce an Auto-MAC matrix close or equal to the identity matrix. Off diagonal elements with considerable magnitudes indicate some trouble with the mode shapes (their orthogonality is compromised) and test and identification methods should be reevaluated.
%% Modal Assurance Criterion (Auto_MAC) for fitted dynamic system
AUTO_MAC_ms_sys_DOF_Hits=AutoMAC(ms_sys_DOF_Hits); %Press CTRL+Enter to re run section.
%

Mode shapes comparison

A series of simulation results were previously computed and will be used for cross comparison. This section compares mode shapes obtained from tests where the impact hammer hits on a given DOF.
For instance mode shape 1 (vibrating at approximately 31.4Hz) should be the same no matter what DOF is hit. hence the MAC value of said mode identified from hits on DOF 3 and 6 should be 1. Additionally, since the true mode shapes of the plate must be orthogonal to each other, comparisons between mode shape 1 of one test and any other mode shape different from 1 should yield a MAC value of 0.
%% Mode shapes comparator %Press CTRL+ENTER keys to re-animate the sequence
DOF_Hammer_1=3;mode_1=1;
DOF_Hammer_2=14;mode_2=1;
fn1=load(strcat('fn_sys_DOF',num2str(DOF_Hammer_1),'_SS50_10Iter.mat'));
dr1=load(strcat('dr_sys_DOF',num2str(DOF_Hammer_1),'_SS50_10Iter.mat'));
ms1=load(strcat('ms_sys_DOF',num2str(DOF_Hammer_1),'_SS50_10Iter.mat'));
dr1=dr1.dr_sys_DOF_Hits;fn1=fn1.fn_sys_DOF_Hits;ms1=ms1.ms_sys_DOF_Hits;
fn2=load(strcat('fn_sys_DOF',num2str(DOF_Hammer_2),'_SS50_10Iter.mat'));
dr2=load(strcat('dr_sys_DOF',num2str(DOF_Hammer_2),'_SS50_10Iter.mat'));
ms2=load(strcat('ms_sys_DOF',num2str(DOF_Hammer_2),'_SS50_10Iter.mat'));
dr2=dr2.dr_sys_DOF_Hits;fn2=fn2.fn_sys_DOF_Hits;ms2=ms2.ms_sys_DOF_Hits;
CAD=VioloncellotopplateCADcoordinates;
Title=strcat('Mode shapes comparisson, system identified with hits on DOF ',num2str(DOF_Hammer_1),'(left) vs system identified with hits on DOF ',num2str(DOF_Hammer_2));
CompareTwoModes(mode_1,fn1,dr1,ms1,mode_2,fn2,dr2,ms2,AccelPos,AccelDir,CAD,Title)
%

References

[1] https://www.mathworks.com/help/ident/examples/modal-analysis-of-a-flexible-flying-wing-aircraft.html#FlexibleWingModesExample-6
[2] Theoretical and experimental modal analysis, ISBN 0 86380 208 7
[3] The physics of the Violin, ISBN 978-0-262-52707-1

Acknowledgements

The authors would like to acknowledge the help of the following people:
To all of them, thank you.
Functions
function Hammer_EMA_test_visualization(Test_number,Sensors_output,sensors_units,hammer_input,hammer_units,input_DOF,T_sampling,post)
%Inputs:
%Test_number -> 1,2,3..Just for graph titles, indicates which test of all
% the ones performed is used for analysis.
%
%Sensors_output -> mxn array containing response signals from n sensors,
% with m data points each.
%
%sensors_units -> 'string', output sensor units i.e. 'm/s^2' or 'g'
%
%hammer_input -> mx1 array containing modal hammer input signal with m data
% points
%input_DOF -> DOF direction of the hammer impact i.e. 1,2,3,... or
% -1,-2,-3... (opposite to sensor positive direction)
%
%T_sampling -> Sampling time in seconds
%Creating time vector
t = (0:length(hammer_input)-1)' * T_sampling;
%Creating figure with two subplots, indicating the test number
fig =figure('Name',strcat('Measured Data Test#',num2str(Test_number)));
set(fig, 'Position', [0 0 1000 500])
subplot(211) %two vertically stacked subplots
plot(t,Sensors_output) %Ploting sensor outputs vs time
title(strcat('Test #',num2str(Test_number),' Output responses, impact hammer on DOF#',num2str(input_DOF)))
legend('Outputs acceleration signals')
ylabel(strcat('Outputs (',sensors_units,')'))
%Looking for the time of maximum peak in the hammer impact
[~,indexOfFirstMax] = max(hammer_input); % Find index of max hammer force value.
maxforcetime = t(indexOfFirstMax);
% Pretrigger
if maxforcetime-post/10>=0
time_pre=maxforcetime-post/10;
else
time_pre=0;
end
xlim([time_pre time_pre+post])
subplot(212)
plot(t,hammer_input)
title(strcat('Test #',num2str(Test_number),' input force, impact hammer on DOF#',num2str(input_DOF)))
legend('hammer impact force')
ylabel(strcat('Input DOF #',num2str(input_DOF),' (',hammer_units,')'))
xlabel('Time (seconds)')
xlim([time_pre time_pre+post])
end
%
%% Plotting FRF estimation for DOF 3 hits
function plot_FRF_Bode(FRF,Output_DOFs,freq_range)
% Inputs;
%
% FRF --> FRF object as created by Matlab's Etfe function
%
% Output_Dofs --> Vector containing the DOFs to plot, i.e. [1 2 15] for
% plotting DOFs 1,2 and 15.
%
% freq_range --> Vector containing the frequency range of the bode plot
opt = bodeoptions; % plot options
opt.FreqUnits = 'Hz'; % show frequency in Hz
opt.PhaseVisible = 'off'; % do not show phase
fig=figure;
set(fig, 'Position', [0 0 1000 500]);
if isempty(Output_DOFs)
bodeplot(FRF, opt); % plot frequency response
else
bodeplot(FRF(Output_DOFs,:), opt) %outputs for plotting
end
if isempty(freq_range)
else
xlim(freq_range)
end
grid on
end
%
%%
function CompareTwoModes(Modenum1,fn1,dr1,Modeshapes1,Modenum2,fn2,dr2,Modeshapes2,AccelPos,AccelDir,CAD,Title)
%Mode shape 1
wn1 = fn1(Modenum1)*2*pi; % Mode frequency in rad/sec
T1 = 1/fn1(Modenum1); % Period of modal oscillation
%Mode shape 2
wn2 = fn2(Modenum2)*2*pi; % Mode frequency in rad/sec
T2 = 1/fn2(Modenum2); % Period of modal oscillation
Np = 4; % Number of periods to simulate
tmax1 = Np*T1; % Simulate Np periods mode 1
tmax2 = Np*T2; % Simulate Np periods mode 1
Nt = 60; % Number of time steps for animation
t1 = linspace(0,tmax1,Nt);
t2 = linspace(0,tmax2,Nt);
ThisMode1 = 0.1*Modeshapes1(:,Modenum1)/max(abs(Modeshapes1(:,Modenum1))); % normalized for plotting
ThisMode2 = 0.1*Modeshapes2(:,Modenum2)/max(abs(Modeshapes2(:,Modenum2))); % normalized for plotting
ThisMode1_unitary=Modeshapes1(:,Modenum1)/norm(Modeshapes1(:,Modenum1));
ThisMode2_unitary=Modeshapes2(:,Modenum2)/norm(Modeshapes2(:,Modenum2));
%Modal Assurance Criterion
MAC=norm(ThisMode1_unitary'*ThisMode2_unitary)^2;
for i=1:length(AccelDir(:,1))
AccelDir(i,:)=AccelDir(i,:)/norm(AccelDir(i,:));
end
fig = figure;
set(fig, 'Position', [0 0 1000 500]);
for k = 1:Nt
Rot1_1 = cos(wn1*t1(k));
Rot2_1 = sin(wn1*t1(k));
%ModeShape1
x_t1 = (AccelDir(:,1)*(real(ThisMode1)*Rot1_1 - imag(ThisMode1)*Rot2_1)')';
x_t1 = AccelPos(:,1)+x_t1(:,1);
y_t1 = (AccelDir(:,2)*(real(ThisMode1(:))*Rot1_1 - imag(ThisMode1(:))*Rot2_1)')';
y_t1 = AccelPos(:,2)+y_t1(:,1);
z_t1 = (AccelDir(:,3)*(real(ThisMode1(:))*Rot1_1 - imag(ThisMode1(:))*Rot2_1)')';
z_t1 = AccelPos(:,3)+z_t1(:,1);
Rot1_2 = cos(wn2*t2(k));
Rot2_2 = sin(wn2*t2(k));
%ModeShape2
x_t2 = (AccelDir(:,1)*(real(ThisMode2)*Rot1_2 - imag(ThisMode2)*Rot2_2)')';
x_t2 = AccelPos(:,1)+x_t2(:,1);
y_t2 = (AccelDir(:,2)*(real(ThisMode2(:))*Rot1_2 - imag(ThisMode2(:))*Rot2_2)')';
y_t2 = AccelPos(:,2)+y_t2(:,1);
z_t2 = (AccelDir(:,3)*(real(ThisMode2(:))*Rot1_2 - imag(ThisMode2(:))*Rot2_2)')';
z_t2 = AccelPos(:,3)+z_t2(:,1);
%Mode shapes extrapolation
F1=scatteredInterpolant(x_t1,y_t1,z_t1);
Vq1=F1(CAD(:,1),CAD(:,2));
F2=scatteredInterpolant(x_t2,y_t2,z_t2);
Vq2=F2(CAD(:,1),CAD(:,2));
pause(0.075)
% Animated by updating mode shape coordinates.
%trimesh(Triangles,AccelPos(:,1),AccelPos(:,2),AccelPos(:,3),'LineStyle','--','FaceAlpha',0.2,'EdgeColor','k');
%Mode Shape 1
subplot(1,2,1)
plot3(CAD(:,1),CAD(:,2),CAD(:,3),'.')
view(60,45)
zlim([-0.5 0.5])
xlim([-0.5 0.5])
ylim([-.1 .9])
hold on
plot3(CAD(:,1),CAD(:,2),Vq1,'o')
title(sprintf('Mode %d. F_N = %.2f Hz damping ratio %.3f MAC= %.3f',Modenum1,fn1(Modenum1),dr1(Modenum1),MAC));
%trimesh(Triangles,x_t,y_t,z_t,'FaceColor','b','LineStyle','-','FaceAlpha',0.1)
view(60,45)
zlim([-0.5 0.5])
xlim([-0.5 0.5])
ylim([-.1 .9])
hold off
%Mode Shape 2
subplot(1,2,2)
plot3(CAD(:,1),CAD(:,2),CAD(:,3),'.')
view(60,45)
zlim([-0.5 0.5])
xlim([-0.5 0.5])
ylim([-.1 .9])
hold on
plot3(CAD(:,1),CAD(:,2),Vq2,'o')
title(sprintf('Mode %d. F_N = %.2f Hz damping ratio %.3f MAC= %.3f', Modenum2,fn2(Modenum2),dr2(Modenum2),MAC));
%trimesh(Triangles,x_t,y_t,z_t,'FaceColor','b','LineStyle','-','FaceAlpha',0.1)
view(60,45)
zlim([-0.5 0.5])
xlim([-0.5 0.5])
ylim([-.1 .9])
hold off
sgtitle(Title)
end
end
function AnimateOneMode(ModeNum, fn,dr, ModeShapes, AccelPos, AccelDir,CAD)
%Delaunay triangularization for accelerometer points
%Triangles=delaunay(AccelPos(:,1),AccelPos(:,2));
%trimesh(T,AccelPos(:,1),AccelPos(:,2),AccelPos(:,3)+0.5);
wn = fn(ModeNum)*2*pi; % Mode frequency in rad/sec
T = 1/fn(ModeNum); % Period of modal oscillation
Np = 4; % Number of periods to simulate
tmax = Np*T; % Simulate Np periods
Nt = 60; % Number of time steps for animation
t = linspace(0,tmax,Nt);
ThisMode = 0.15*ModeShapes(:,ModeNum)/max(abs(ModeShapes(:,ModeNum))); % normalized for plotting
for i=1:length(AccelDir(:,1))
AccelDir(i,:)=AccelDir(i,:)/norm(AccelDir(i,:));
end
fig=figure;
set(fig, 'Position', [0 0 1000 500]);
for k = 1:Nt
Rot1 = cos(wn*t(k));
Rot2 = sin(wn*t(k));
x_t = (AccelDir(:,1)*(real(ThisMode)*Rot1 - imag(ThisMode)*Rot2)')';
x_t = AccelPos(:,1)+x_t(:,1);
y_t = (AccelDir(:,2)*(real(ThisMode(:))*Rot1 - imag(ThisMode(:))*Rot2)')';
y_t = AccelPos(:,2)+y_t(:,1);
z_t = (AccelDir(:,3)*(real(ThisMode(:))*Rot1 - imag(ThisMode(:))*Rot2)')';
z_t = AccelPos(:,3)+z_t(:,1);
F=scatteredInterpolant(x_t,y_t,z_t);
Vq=F(CAD(:,1),CAD(:,2));
pause(0.075)
% Animated by updating mode shape coordinates.
%trimesh(Triangles,AccelPos(:,1),AccelPos(:,2),AccelPos(:,3),'LineStyle','--','FaceAlpha',0.2,'EdgeColor','k');
plot3(CAD(:,1),CAD(:,2),CAD(:,3),'.')
view(20,60)
zlim([-0.5 0.5])
xlim([-0.5 0.5])
ylim([-.1 .9])
hold on
plot3(CAD(:,1),CAD(:,2),Vq,'o')
title(sprintf('Mode %d. FN = %.1f Hz damping ratio %.3f', ModeNum,fn(ModeNum),dr(ModeNum)));
%trimesh(Triangles,x_t,y_t,z_t,'FaceColor','b','LineStyle','-','FaceAlpha',0.1)
view(20,60)
zlim([-0.5 0.5])
xlim([-0.5 0.5])
ylim([-.1 .9])
hold off
end
end
function [Matrix]=AutoMAC(ms_array)
ms_array_normalized=zeros(size(ms_array));
for i=1:length(ms_array(1,:))
ms_array_normalized(:,i)=ms_array(:,i)/norm(ms_array(:,i));
end
Matrix=zeros(size(ms_array'*ms_array));
for i=1:length(ms_array(1,:))
for j=1:length(ms_array(1,:))
Matrix(i,j)= norm(ms_array_normalized(:,j)'*ms_array_normalized(:,i))^2;
end
end
bar3c(Matrix)
colorbar
end
function hh = bar3c( varargin )
%BAR3C Extension of bar3, which sets the bar color corresponding to its
%height.
%
% Extra parameter: 'MaxColorVal', maxclr
% This will make the color/height of the bar absolute against this maximum
% value. Otherwise, the colors will be relative against the maximum zdata
% value of all bars on the graph.
%
[abscolor, idxabsc]=getarg('MaxColorVal',varargin{:});
if idxabsc
varargin(idxabsc+(0:1))=[];
end
h = bar3(varargin{:});
for ii = 1:numel(h)
zdata = get(h(ii),'Zdata');
cdata = makecdata(zdata(2:6:end,2),abscolor);
set(h(ii),'Cdata',cdata, 'facecolor','flat');
end
if nargout>0,
hh = h;
end
end
function [val, idx] = getarg(strname,varargin)
idx = 0;
val = NaN;
for jj=1:nargin-2
if strcmpi(varargin{jj},strname)
idx = jj;
val = varargin{jj+1};
return;
end
end
end
function cdata = makecdata(clrs,maxclr)
cdata = NaN(6*numel(clrs),4);
for ii=1:numel(clrs)
cdata((ii-1)*6+(1:6),:)=makesingle_cdata(clrs(ii));
end
if nargin>=2
% it doesn't matter that we put an extra value on cdata(1,1)
% this vertex is NaN (serves as a separator
cdata(1,1)=maxclr;
end
end
function scdata = makesingle_cdata(clr)
scdata = NaN(6,4);
scdata(sub2ind(size(scdata),[3,2,2,1,2,4],[2,1,2,2,3,2]))=clr;
end